data <- readRDS("very_low_birthweight.RDS")
Загрузите датасет very_low_birthweight.RDS (лежит в папке домашнего
задания). Это данные о 671 младенце с очень низкой массой тела (<1600
грамм), собранные в Duke University Medical Center доктором Майклом О’Ши
c 1981 по 1987 г.
Описание переменных см. здесь. Переменными исхода являются колонки
‘dead’, а также время от рождения до смерти или выписки (выводятся из
‘birth’ и ‘exit’ 7 пациентов были выписаны до рождения). Сделайте копию
датасета, в которой удалите колонки с количеством пропусков больше 100,
а затем удалите все строки с пропусками.
data %>%
head()
## birth exit hospstay lowph pltct race bwt gest inout twn lol
## 1 81.511 81.604 34 NA 100 white 1250 35 born at Duke 0 NA
## 2 81.514 81.539 9 7.250000 244 white 1370 32 born at Duke 0 NA
## 3 81.552 81.552 -2 7.059998 114 black 620 23 born at Duke 0 NA
## 4 81.558 81.667 40 7.250000 182 black 1480 32 born at Duke 0 NA
## 5 81.593 81.599 2 6.969997 54 black 925 28 born at Duke 0 NA
## 6 81.602 81.771 62 7.189999 NA white 940 28 born at Duke 0 NA
## magsulf meth toc delivery apg1 vent pneumo pda cld pvh ivh ipe
## 1 NA 0 0 abdominal 8 0 0 0 0 <NA> <NA> <NA>
## 2 NA 1 0 abdominal 7 0 0 0 0 <NA> <NA> <NA>
## 3 NA 0 1 vaginal 1 1 0 0 NA <NA> <NA> <NA>
## 4 NA 1 0 vaginal 8 0 0 0 0 <NA> <NA> <NA>
## 5 NA 0 0 abdominal 5 1 1 0 0 definite definite <NA>
## 6 NA 1 0 abdominal 8 1 0 0 0 absent absent absent
## year sex dead
## 1 81.51196 female 0
## 2 81.51471 female 0
## 3 81.55304 female 1
## 4 81.55847 male 0
## 5 81.59406 female 1
## 6 81.60229 female 0
cols_to_select <- data.frame(NA_count = colSums(is.na(data))) %>% filter(NA_count <= 100) %>% row.names()
data %>%
select(all_of(cols_to_select)) %>%
na.omit() %>%
write_csv("filtered_dataset.csv")
Постройте графики плотности распределения для числовых переменных. Удалите выбросы, если таковые имеются. Преобразуйте категориальные переменные в факторы. Для любых двух числовых переменных раскрасьте график по переменной ‘inout’.
data <- data %>%
mutate(across(where(is.character), as.numeric),
across(where(is.factor), as.factor),
across(where(~ is.numeric(.) & n_distinct(., na.rm = TRUE) == 2), as.factor))
data %>%
select_if(is.numeric) %>%
pivot_longer(cols = everything(), names_to = "variable", values_to = "value") %>%
ggplot(aes(x = value)) +
geom_histogram(bins = 10, fill = "skyblue", color = "black") +
facet_wrap(~variable, scales = "free_x") +
theme_minimal() +
theme(axis.text.x = element_text(size=7))
## Warning: Removed 657 rows containing non-finite outside the scale range
## (`stat_bin()`).
dfSummary(data)
## Data Frame Summary
## data
## Dimensions: 671 x 26
## Duplicates: 0
##
## --------------------------------------------------------------------------------------------------------------
## No Variable Stats / Values Freqs (% of Valid) Graph Valid Missing
## ---- ----------- ---------------------------- --------------------- --------------------- ---------- ---------
## 1 birth Mean (sd) : 84.8 (1.6) 542 distinct values : : 650 21
## [numeric] min < med < max: : : : : . : (96.9%) (3.1%)
## 81.5 < 84.9 < 87.5 : : : : : : : : :
## IQR (CV) : 2.6 (0) : : : : : : : : : :
## : : : : : : : : : :
##
## 2 exit Mean (sd) : 84.8 (1.8) 539 distinct values : 640 31
## [numeric] min < med < max: : . (95.4%) (4.6%)
## 68.5 < 85 < 96.9 : :
## IQR (CV) : 2.6 (0) : :
## : : :
##
## 3 hospstay Mean (sd) : 40.4 (304.8) 153 distinct values : 640 31
## [integer] min < med < max: : (95.4%) (4.6%)
## -6574 < 37 < 3668 :
## IQR (CV) : 46 (7.6) :
## :
##
## 4 lowph Mean (sd) : 7.2 (0.1) 74 distinct values : . 609 62
## [numeric] min < med < max: : : (90.8%) (9.2%)
## 6.5 < 7.2 < 7.5 . : :
## IQR (CV) : 0.2 (0) . : : : :
## . : : : : : .
##
## 5 pltct Mean (sd) : 201.6 (80.5) 266 distinct values : 601 70
## [integer] min < med < max: . : . (89.6%) (10.4%)
## 16 < 202 < 571 . : : :
## IQR (CV) : 109 (0.4) : : : : .
## : : : : : : .
##
## 6 race 1. black 369 (57.1%) IIIIIIIIIII 646 25
## [factor] 2. native American 16 ( 2.5%) (96.3%) (3.7%)
## 3. oriental 4 ( 0.6%)
## 4. white 257 (39.8%) IIIIIII
##
## 7 bwt Mean (sd) : 1093.9 (265.2) 142 distinct values . : 669 2
## [integer] min < med < max: . : : : : (99.7%) (0.3%)
## 400 < 1120 < 1580 : : : : : :
## IQR (CV) : 410 (0.2) . : : : : : : : :
## . : : : : : : : : :
##
## 8 gest Mean (sd) : 28.9 (2.5) 24 distinct values : : 667 4
## [numeric] min < med < max: : : : (99.4%) (0.6%)
## 22 < 29 < 40 . : : :
## IQR (CV) : 4 (0.1) : : : :
## . : : : : :
##
## 9 inout 1. born at Duke 547 (81.9%) IIIIIIIIIIIIIIII 668 3
## [factor] 2. transported 121 (18.1%) III (99.6%) (0.4%)
##
## 10 twn 1. 0 516 (79.3%) IIIIIIIIIIIIIII 651 20
## [factor] 2. 1 135 (20.7%) IIII (97.0%) (3.0%)
##
## 11 lol Mean (sd) : 8.4 (19.3) 40 distinct values : 290 381
## [integer] min < med < max: : (43.2%) (56.8%)
## 0 < 3.5 < 192 :
## IQR (CV) : 9 (2.3) :
## : .
##
## 12 magsulf 1. 0 367 (86.6%) IIIIIIIIIIIIIIIII 424 247
## [factor] 2. 1 57 (13.4%) II (63.2%) (36.8%)
##
## 13 meth 1. 0 318 (56.3%) IIIIIIIIIII 565 106
## [factor] 2. 1 247 (43.7%) IIIIIIII (84.2%) (15.8%)
##
## 14 toc 1. 0 438 (77.5%) IIIIIIIIIIIIIII 565 106
## [factor] 2. 1 127 (22.5%) IIII (84.2%) (15.8%)
##
## 15 delivery 1. abdominal 314 (48.4%) IIIIIIIII 649 22
## [factor] 2. vaginal 335 (51.6%) IIIIIIIIII (96.7%) (3.3%)
##
## 16 apg1 Mean (sd) : 4.9 (2.6) 0 : 5 ( 0.8%) 637 34
## [integer] min < med < max: 1 : 92 (14.4%) II (94.9%) (5.1%)
## 0 < 5 < 9 2 : 69 (10.8%) II
## IQR (CV) : 5 (0.5) 3 : 57 ( 8.9%) I
## 4 : 45 ( 7.1%) I
## 5 : 69 (10.8%) II
## 6 : 80 (12.6%) II
## 7 : 80 (12.6%) II
## 8 : 103 (16.2%) III
## 9 : 37 ( 5.8%) I
##
## 17 vent 1. 0 269 (42.0%) IIIIIIII 641 30
## [factor] 2. 1 372 (58.0%) IIIIIIIIIII (95.5%) (4.5%)
##
## 18 pneumo 1. 0 518 (80.3%) IIIIIIIIIIIIIIII 645 26
## [factor] 2. 1 127 (19.7%) III (96.1%) (3.9%)
##
## 19 pda 1. 0 508 (79.1%) IIIIIIIIIIIIIII 642 29
## [factor] 2. 1 134 (20.9%) IIII (95.7%) (4.3%)
##
## 20 cld 1. 0 442 (73.1%) IIIIIIIIIIIIII 605 66
## [factor] 2. 1 163 (26.9%) IIIII (90.2%) (9.8%)
##
## 21 pvh 1. absent 360 (68.4%) IIIIIIIIIIIII 526 145
## [factor] 2. definite 125 (23.8%) IIII (78.4%) (21.6%)
## 3. possible 41 ( 7.8%) I
##
## 22 ivh 1. absent 442 (83.9%) IIIIIIIIIIIIIIII 527 144
## [factor] 2. definite 75 (14.2%) II (78.5%) (21.5%)
## 3. possible 10 ( 1.9%)
##
## 23 ipe 1. absent 472 (89.6%) IIIIIIIIIIIIIIIII 527 144
## [factor] 2. definite 38 ( 7.2%) I (78.5%) (21.5%)
## 3. possible 17 ( 3.2%)
##
## 24 year Mean (sd) : 84.8 (1.6) 542 distinct values : : 650 21
## [numeric] min < med < max: : : : : . : (96.9%) (3.1%)
## 81.5 < 84.9 < 87.5 : : : : : : : : :
## IQR (CV) : 2.6 (0) : : : : : : : : : :
## : : : : : : : : : :
##
## 25 sex 1. female 320 (49.2%) IIIIIIIII 650 21
## [factor] 2. male 330 (50.8%) IIIIIIIIII (96.9%) (3.1%)
##
## 26 dead 1. 0 527 (78.5%) IIIIIIIIIIIIIII 671 0
## [factor] 2. 1 144 (21.5%) IIII (100.0%) (0.0%)
## --------------------------------------------------------------------------------------------------------------
Что-то странное происходит в hospstay, предполагаю, что младенцы не умеют путешествовать во времени. В этих (exit, pltct, lol) переменных так же есть экстремальные значения, природа которых мне неизвестна, поэтому удалю их по правилу трех IQR.
NB: после удаления аутлаеров по hospstay исчезли аутлаеры у exit
data <- data %>%
mutate(ID = c(1:nrow(data))) %>%
select(ID, everything())
data <- data %>%
filter(hospstay > 0 | is.na(hospstay),
hospstay < 3 * IQR(hospstay, na.rm = T) | is.na(hospstay),
pltct < 3 * IQR(pltct, na.rm = T) | is.na(pltct),
lol < 3 * IQR(lol, na.rm = T) | is.na(lol))
data %>%
dfSummary()
## Data Frame Summary
## Dimensions: 569 x 27
## Duplicates: 0
##
## ---------------------------------------------------------------------------------------------------------------
## No Variable Stats / Values Freqs (% of Valid) Graph Valid Missing
## ---- ----------- ---------------------------- --------------------- ---------------------- ---------- ---------
## 1 ID Mean (sd) : 332.9 (193.8) 569 distinct values . : : . : . . . . . 569 0
## [integer] min < med < max: : : : : : : : : : : (100.0%) (0.0%)
## 1 < 330 < 671 : : : : : : : : : :
## IQR (CV) : 336 (0.6) : : : : : : : : : :
## : : : : : : : : : :
##
## 2 birth Mean (sd) : 84.7 (1.6) 460 distinct values : . 548 21
## [numeric] min < med < max: : : : : . (96.3%) (3.7%)
## 81.5 < 84.9 < 87.5 : : : : : : : : :
## IQR (CV) : 2.5 (0) : : : : : : : : : :
## : : : : : : : : : :
##
## 3 exit Mean (sd) : 84.8 (1.6) 460 distinct values : 538 31
## [numeric] min < med < max: : . : : : (94.6%) (5.4%)
## 81.5 < 84.9 < 87.8 : : : : : : : : .
## IQR (CV) : 2.6 (0) . : : : : : : : : :
## : : : : : : : : : :
##
## 4 hospstay Mean (sd) : 41.6 (30.6) 117 distinct values : 538 31
## [integer] min < med < max: : : : . (94.6%) (5.4%)
## 1 < 36 < 137 : : : : .
## IQR (CV) : 43 (0.7) : : : : : .
## : : : : : : : . . .
##
## 5 lowph Mean (sd) : 7.2 (0.1) 71 distinct values : . 516 53
## [numeric] min < med < max: : : (90.7%) (9.3%)
## 6.5 < 7.2 < 7.5 . : :
## IQR (CV) : 0.2 (0) . : : : :
## . : : : : : .
##
## 6 pltct Mean (sd) : 192.1 (68.1) 222 distinct values : 509 60
## [integer] min < med < max: : : (89.5%) (10.5%)
## 16 < 199 < 324 : : : .
## IQR (CV) : 99 (0.4) . : : : :
## . : : : : : :
##
## 7 race 1. black 303 (55.7%) IIIIIIIIIII 544 25
## [factor] 2. native American 16 ( 2.9%) (95.6%) (4.4%)
## 3. oriental 3 ( 0.6%)
## 4. white 222 (40.8%) IIIIIIII
##
## 8 bwt Mean (sd) : 1110.1 (258.9) 134 distinct values : 567 2
## [integer] min < med < max: . : : : (99.6%) (0.4%)
## 400 < 1150 < 1580 . : : : : :
## IQR (CV) : 400 (0.2) . : : : : : : .
## . : : : : : : : : :
##
## 9 gest Mean (sd) : 29.1 (2.5) 23 distinct values . : 566 3
## [numeric] min < med < max: : : : (99.5%) (0.5%)
## 23 < 29 < 40 : : :
## IQR (CV) : 4 (0.1) : : : :
## . : : : : :
##
## 10 inout 1. born at Duke 462 (81.5%) IIIIIIIIIIIIIIII 567 2
## [factor] 2. transported 105 (18.5%) III (99.6%) (0.4%)
##
## 11 twn 1. 0 432 (78.5%) IIIIIIIIIIIIIII 550 19
## [factor] 2. 1 118 (21.5%) IIII (96.7%) (3.3%)
##
## 12 lol Mean (sd) : 5.2 (6.3) 26 distinct values : 242 327
## [integer] min < med < max: : (42.5%) (57.5%)
## 0 < 3 < 26 :
## IQR (CV) : 8 (1.2) : .
## : : . . .
##
## 13 magsulf 1. 0 306 (85.2%) IIIIIIIIIIIIIIIII 359 210
## [factor] 2. 1 53 (14.8%) II (63.1%) (36.9%)
##
## 14 meth 1. 0 263 (55.4%) IIIIIIIIIII 475 94
## [factor] 2. 1 212 (44.6%) IIIIIIII (83.5%) (16.5%)
##
## 15 toc 1. 0 372 (78.5%) IIIIIIIIIIIIIII 474 95
## [factor] 2. 1 102 (21.5%) IIII (83.3%) (16.7%)
##
## 16 delivery 1. abdominal 274 (50.1%) IIIIIIIIII 547 22
## [factor] 2. vaginal 273 (49.9%) IIIIIIIII (96.1%) (3.9%)
##
## 17 apg1 Mean (sd) : 5 (2.6) 0 : 3 ( 0.6%) 536 33
## [integer] min < med < max: 1 : 74 (13.8%) II (94.2%) (5.8%)
## 0 < 5 < 9 2 : 53 ( 9.9%) I
## IQR (CV) : 4 (0.5) 3 : 46 ( 8.6%) I
## 4 : 35 ( 6.5%) I
## 5 : 61 (11.4%) II
## 6 : 68 (12.7%) II
## 7 : 71 (13.2%) II
## 8 : 96 (17.9%) III
## 9 : 29 ( 5.4%) I
##
## 18 vent 1. 0 236 (43.5%) IIIIIIII 543 26
## [factor] 2. 1 307 (56.5%) IIIIIIIIIII (95.4%) (4.6%)
##
## 19 pneumo 1. 0 441 (81.1%) IIIIIIIIIIIIIIII 544 25
## [factor] 2. 1 103 (18.9%) III (95.6%) (4.4%)
##
## 20 pda 1. 0 441 (81.1%) IIIIIIIIIIIIIIII 544 25
## [factor] 2. 1 103 (18.9%) III (95.6%) (4.4%)
##
## 21 cld 1. 0 393 (76.2%) IIIIIIIIIIIIIII 516 53
## [factor] 2. 1 123 (23.8%) IIII (90.7%) (9.3%)
##
## 22 pvh 1. absent 307 (69.3%) IIIIIIIIIIIII 443 126
## [factor] 2. definite 100 (22.6%) IIII (77.9%) (22.1%)
## 3. possible 36 ( 8.1%) I
##
## 23 ivh 1. absent 378 (85.3%) IIIIIIIIIIIIIIIII 443 126
## [factor] 2. definite 57 (12.9%) II (77.9%) (22.1%)
## 3. possible 8 ( 1.8%)
##
## 24 ipe 1. absent 399 (90.1%) IIIIIIIIIIIIIIIIII 443 126
## [factor] 2. definite 30 ( 6.8%) I (77.9%) (22.1%)
## 3. possible 14 ( 3.2%)
##
## 25 year Mean (sd) : 84.7 (1.6) 460 distinct values : . 548 21
## [numeric] min < med < max: : : : : . (96.3%) (3.7%)
## 81.5 < 84.9 < 87.5 : . : : : : : : :
## IQR (CV) : 2.5 (0) . : : : : : : : : :
## : : : : : : : : : :
##
## 26 sex 1. female 277 (50.5%) IIIIIIIIII 548 21
## [factor] 2. male 271 (49.5%) IIIIIIIII (96.3%) (3.7%)
##
## 27 dead 1. 0 462 (81.2%) IIIIIIIIIIIIIIII 569 0
## [factor] 2. 1 107 (18.8%) III (100.0%) (0.0%)
## ---------------------------------------------------------------------------------------------------------------
data %>%
select_if(is.numeric) %>%
select(-ID) %>%
pivot_longer(cols = everything(), names_to = "variable", values_to = "value") %>%
ggplot(aes(x = value)) +
geom_histogram(bins = 10, fill = "skyblue", color = "black") +
facet_wrap(~variable, scales = "free_x") +
theme_minimal() +
theme(axis.text.x = element_text(size=7))
## Warning: Removed 582 rows containing non-finite outside the scale range
## (`stat_bin()`).
Теперь выглядит приятнее :)
Проведите тест на сравнение значений колонки ‘lowph’ между группами в переменной inout. Вид статистического теста определите самостоятельно. Визуализируйте результат через библиотеку ‘rstatix’. Как бы вы интерпретировали результат, если бы знали, что более низкое значение lowph ассоциировано с более низкой выживаемостью?
data %>% group_by(inout) %>% select(inout, lowph) %>% dfSummary()
## Data Frame Summary
## Group: inout = born at Duke
## Dimensions: 462 x 2
## Duplicates: 396
##
## -------------------------------------------------------------------------------------------------------
## No Variable Stats / Values Freqs (% of Valid) Graph Valid Missing
## ---- ----------- ----------------------- -------------------- --------------------- --------- ---------
## 2 lowph Mean (sd) : 7.2 (0.1) 65 distinct values : : 420 42
## [numeric] min < med < max: : : (90.9%) (9.1%)
## 6.5 < 7.2 < 7.5 . : : .
## IQR (CV) : 0.2 (0) : : : :
## . : : : : : .
## -------------------------------------------------------------------------------------------------------
##
## Group: inout = transported
## Dimensions: 105 x 2
## Duplicates: 60
##
## ---------------------------------------------------------------------------------------------------
## No Variable Stats / Values Freqs (% of Valid) Graph Valid Missing
## ---- ----------- ----------------------- -------------------- ----------------- --------- ---------
## 2 lowph Mean (sd) : 7.1 (0.2) 44 distinct values : 96 9
## [numeric] min < med < max: : . (91.4%) (8.6%)
## 6.7 < 7.2 < 7.5 : :
## IQR (CV) : 0.2 (0) : : : : .
## . . . : : : : :
## ---------------------------------------------------------------------------------------------------
##
## Group: inout = NA
## Dimensions: 2 x 2
## Duplicates: 1
##
## ----------------------------------------------------------------------------------
## No Variable Stats / Values Freqs (% of Valid) Graph Valid Missing
## ---- ----------- ---------------- -------------------- ------- -------- ----------
## 2 lowph All NA's 0 2
## [numeric] (0.0%) (100.0%)
## ----------------------------------------------------------------------------------
data %>%
t_test(lowph ~ inout, detailed = T)
## Warning: The statistic is based on a difference or ratio; by default, for
## difference-based statistics, the explanatory variable is subtracted in the
## order "born at Duke" - "transported", or divided in the order "born at Duke" /
## "transported" for ratio-based statistics. To specify this order yourself,
## supply `order = c("born at Duke", "transported")`.
## # A tibble: 1 × 7
## statistic t_df p_value alternative estimate lower_ci upper_ci
## <dbl> <dbl> <dbl> <chr> <dbl> <dbl> <dbl>
## 1 5.20 130. 0.000000744 two.sided 0.0881 0.0546 0.122
У транспортированных пациентов значение рН меньше, чем у рожденных в Дюке (MD=0.088, CI 95%: [0.055; 0.122], p < \(10^{-6}\)). Поскольку мы знаем, что более низкое значение рН ассоциировано с более высокой смертностью, можно предположить, что транспортировка пациентов может быть также связана с более высокой смертностью.
Сделайте новый датафрейм, в котором оставьте только континуальные или ранговые данные, кроме ‘birth’, ‘year’ и ‘exit’. Сделайте корреляционный анализ этих данных. Постройте два любых типа графиков для визуализации корреляций.
corr_data <- data %>%
select(where(is.numeric), -c(ID, birth, year, exit)) %>%
psych::corr.test(method = "spearman", adjust="BH")
cor_melted <- reshape2::melt(corr_data$r)
p_melted <- reshape2::melt(corr_data$p)
cor_plot_data <- merge(cor_melted, p_melted, by = c("Var1", "Var2"))
colnames(cor_plot_data) <- c("Var1", "Var2", "Correlation", "Adjusted_P")
ggplot(cor_plot_data, aes(Var1, Var2)) +
geom_tile(aes(fill = Correlation), color = "white") +
scale_fill_gradient2(low = "#433CE1", high = "#E13C59", mid = "grey90",
midpoint = 0, limit = c(-1, 1), name = "Correlation") +
geom_text(aes(label = ifelse(Adjusted_P < 0.05, sprintf("%.2f", Adjusted_P), "")),
color = "black", size = 4) +
theme_minimal() +
theme(axis.text.x = element_text(angle = 45, hjust = 1, size = 12),
axis.text.y = element_text(size = 12)) +
labs(x = "", y = "")
data %>%
select(where(is.numeric), -c(ID, birth, year, exit)) %>%
ggpairs(progress = FALSE) +
theme_bw() +
theme(axis.text = element_text(size=6))
Постройте иерархическую кластеризацию на этом датафрейме.
rownames(data) <- data$ID
res <- data %>%
select(where(is.numeric), -c(ID, birth, year, exit)) %>%
na.omit() %>%
scale() %>%
dist(method = "euclidean") %>%
hclust(method = "ward.D2")
res %>%
fviz_dend(cex = 0.2,
k = 3,
k_colors = "jco")
## Warning: The `<scale>` argument of `guides()` cannot be `FALSE`. Use "none" instead as
## of ggplot2 3.3.4.
## ℹ The deprecated feature was likely used in the factoextra package.
## Please report the issue at <https://github.com/kassambara/factoextra/issues>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
Сделайте одновременный график heatmap и иерархической кластеризации. Интерпретируйте результат.
data %>%
select(where(is.numeric), -c(ID, birth, year, exit)) %>%
na.omit() %>%
scale() %>%
pheatmap::pheatmap(color=colorRampPalette(c("#433CE1", "grey90", "#E13C59"))(50),
fontsize_row = 0.00001, fontsize_col = 7, angle_col = 0)
data %>%
select(where(is.numeric), -c(ID, birth, year, exit)) %>%
na.omit() %>%
scale() %>%
pheatmap::pheatmap(color=colorRampPalette(c("#433CE1", "grey90", "#E13C59"))(50),
fontsize_row = 0.00001, fontsize_col = 7, angle_col = 0, kmeans_k=4)
Визульно, можно выделить на хитмапе 1 4 кластера, чтобы их проинтерпретировать, можно усреднить значения в кластере, как на хитмапе 2. Первые 2 кластера довольно сильно похожи друг на друга, за исключением того, что пациенты из 1 кластера имели более выскойи уровень тромбоцитов и скор APG в отличие от пациентов из кластера 1. Пациенты из кластера 3 в отличие от остальных имеют наиболее длительный срок нахождения в больнице. Отличительными чертами пациентов из кластера 4 является самый низкий уровень рН, масса тела при рождении, а также минимальные сроки беременности.
Проведите PCA анализ на этих данных. Проинтерпретируйте результат. Нужно ли применять шкалирование для этих данных перед проведением PCA?
Постройте biplot график для PCA. Раскрасьте его по значению колонки ‘dead’.
pca.data <- data %>%
select(where(is.numeric), -c(ID, birth, year, exit)) %>%
na.omit() %>%
scale() %>%
prcomp()
fviz_eig(pca.data,
addlabels = T,
ylim = c(0, 50))
В нашем случае необходимо проводить шкалирование данных, так как все переменные имеют разные единицы измерения и разный размах. Первые две компоненты довольно хорошо объясняют данные, объясняя около 55% вариабельности.
pca.data$x %>%
ggplot() +
geom_point(aes(x = PC1, y = PC2)) +
theme_bw()
Вклад каждой переменной в компоненту:
pcs <- reshape2::melt(pca.data$rotation)
ggplot(pcs, aes(Var1, Var2)) +
geom_tile(aes(fill = value), color = "white") +
scale_fill_gradient2(low = "#433CE1", high = "#E13C59", mid = "grey90",
midpoint = 0, limit = c(-1, 1)) +
geom_text(aes(label = sprintf("%.2f", value)),
color = "black", size = 4) +
theme_minimal() +
theme(axis.text.x = element_text(angle = 45, hjust = 1, size = 12),
axis.text.y = element_text(size = 12)) +
labs(x = "", y = "")
В первую компоненту основной положительный вклад имеет время,
проведенное в госпитале, и негативная масса тела при рождении и время
беременности (что похоже на последний кластер в анализе выше).
Во вторую компоненту основной вклад вносит продолжительность родов. В третий уровень тромбоцитов и APG скор в большей степени и в меньшей степени продолжительность нахождения в больнице.
deaths <- data %>% filter(ID %in% rownames(pca.data$x)) %>% select(dead)
levels(deaths$dead) <- c("Alive", "Dead")
ggbiplot(pca.data, choices=1:2,
scale=1, alpha = 0.5, groups = as.factor(deaths$dead)) +
theme_bw()
Переведите последний график в ‘plotly’. При наведении на точку нужно, чтобы отображалось id пациента.
explained_variance <- pca.data$sdev
explained_variance_pct <- round(100 * explained_variance ** 2 / sum(explained_variance ** 2), 2)
comp <- pca.data$rotation
loadings <- comp
for (i in seq(explained_variance)){
loadings[,i] <- comp[,i] * explained_variance[i]
}
features <- rownames(pca.data$rotation)
fig <- plot_ly(
data = as.data.frame(pca.data$x), x = ~PC1, y = ~PC2,
text = ~paste("ID: ", rownames(pca.data$x)),
color = ~as.factor(deaths$dead),
colors = c('#636EFA','#EF553B'),
type = 'scatter',
mode = 'markers') %>%
layout(
legend=list(title=list(text='color')),
plot_bgcolor = "#e5ecf6",
xaxis = list(
title = sprintf("standardized PC1 (%.2f%% explained var.)", explained_variance_pct[1])),
yaxis = list(
title = sprintf("standardized PC2 (%.2f%% explained var.)", explained_variance_pct[2])))
for (i in seq(length(features))){
fig <- fig %>%
add_segments(x = 0,
xend = loadings[i, 1],
y = 0,
yend = loadings[i, 2],
line = list(color = 'black'),
inherit = FALSE,
showlegend = FALSE) %>%
add_annotations(x=loadings[i, 1],
y=loadings[i, 2],
ax = 0,
ay = 0,
text = features[i],
xanchor = 'center',
yanchor= 'bottom')
}
fig
Дайте содержательную интерпретацию PCA анализу. Почему использовать колонку ‘dead’ для выводов об ассоциации с выживаемостью некорректно?
To be done
Приведите ваши данные к размерности в две колонки через UMAP. Сравните результаты отображения точек между алгоритмами PCA и UMAP
umap_comps <- data %>%
select(where(is.numeric), -c(ID, birth, year, exit)) %>%
na.omit() %>%
recipe(~.) %>%
step_normalize(all_predictors()) %>%
step_umap(all_predictors()) %>%
prep() %>%
juice()
combined_data <- cbind(umap_comps, pca.data$x[,1:2], deaths)
umap_plot <- combined_data %>%
ggplot(aes(UMAP1, UMAP2)) +
geom_point(aes(color = as.character(dead)),
alpha = 0.7, size = 2) +
labs(color = NULL) +
theme_bw()
pca_plot <- combined_data %>%
ggplot(aes(PC1, PC2)) +
geom_point(aes(color = as.character(dead)),
alpha = 0.7, size = 2) +
labs(color = NULL) +
theme_bw()
ggarrange(
pca_plot, umap_plot,
nrow = 1, ncol = 2,
widths = c(1, 1)
)
В целом, результаты получились похожими, хотя умершие субъекты лучше
“кластеризуются” в UMAP.